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Abstract. We present results from fully nonlinear simulations of unequal mass binary 
black holes plunging from close separations well inside the innermost stable circular 
orbit with mass ratios q = Mi/M 2 = {1,0.85,0.78,0.55,0.32}, or equivalently, with 
reduced mass parameters rj = M 1 M 2 /{M 1 + M 2 ) 2 = {0.25,0.248,0.246,0.229,0.183}. 
For each case, the initial binary orbital parameters are chosen from the Cook- 
Baumgarte equal-mass ISCO configuration. We show waveforms of the dominant 
£ = 2,3 modes and compute estimates of energy and angular momentum radiated. 
For the plunges from the close separations considered, we measure kick velocities 
from gravitational radiation recoil in the range 25-82 km/s. Due to the initial close 
separations our kick velocity estimates should be understood as a lower bound. 
The close configurations considered are also likely to contain significant eccentricities 
influencing the recoil velocity. 



PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 95.30.Sf, 97.60.Lf 
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1. Introduction 

The Laser Interferometer Space Antenna (LISA) will offer a unique look at merging 
supermassive black holes (BHs) [3J, H|. When galaxies collide, the BHs at the center 
of each galaxy [5l [6] will merge and radiate gravitational waves. In a generic merger 
situation, the colliding BHs will have different masses and spins, thus they will radiate 
gravitational waves anisotropically. This anisotropic radiation will carry both a net 
angular and linear momentum [7] . The linear momentum emitted implies a recoil of, or 
kick to, the BH product of the merger. A large enough recoil may cause the BH resulting 
from the merger to be strongly offset from the center of the galaxy or potentially even 
kicked out of small galaxies [8j [9] . Such a scenario would have a significant impact on 
the standard picture of merger tree history of galaxies. 

A complete investigation of the gravitational recoil from unequal mass BH mergers 
must include the contributions to the kick from the inspiral, merger and ringdown stages 
of the life of the binary. Given the tremendous progress that numerical relativity has 
recently made in simulating binary black hole (BBH) systems [TUJ [HI O [131 HH [15] . 
fully numerical investigations can now be used to investigate gravitational recoil. Recent 
numerical studies [36], [18] as well as post-Newtonian estimates [T6l [TT] suggest that most 
of the kick gets accumulated during the plunge. Motivated by these observations, 
we carried out a series of unequal mass BBHs simulations, using initial data that 
corresponds to BHs plunging from innermost stable circular orbit (ISCO), i.e. initial 
separations of 2.34M. In particular, we focus on the dominant I = 2, 3 modes and 
study how the energy radiated in these modes changes as the mass ratio is varied. 
Additionally, we calculate the energy and angular momentum radiated, and from the 
gravitational radiation we estimate kick velocities in the range 25-82 km/s acquired by 
the final BH. 

Soon after the completion of our study, Baker et al. [36J studied a much further 
separated system and Gonzales et al. [18] published fully relativistic kick estimates from 
BHs with initial separation of 7M for a broad range of mass ratios, from q = M1/M2 = 1 
to q = 0.25, or equivalently 77 = q/(l+q) 2 from 0.25 to 0.16. They estimated a maximum 
kick of 175.70 kms" 1 for a mass ratio of 77 = 0.195. Using close-limit approximation 
techniques Sopuerta et al. investigated gravitational wave recoil and also the effect of 
small eccentricities [U [2] . 

Initial Data. The initial data sets are constructed via the puncture method using a 
spectral code [191 120] . The essence of this method is to solve the Hamiltonian constraint 
for the conformal factor <fi. The initial three- metric is conformally flat, maximally sliced, 
and the extrinsic curvature is given by the Bowen-York solution to the momentum 
constraint. The conformal factor <fr is used to set the initial lapse as a = (f)~ 2 [21] 113] . 
while the initial shift is f3 l = 0. 

For the equal mass setup, we evolve the so-called QC-0 initial data set [22]. This 
is intended to represent a quasi-circular configuration of inspiraling puncture BBHs at 
the ISCO. QC-0 data have been used as the starting point by other studies [TTl [T3] [Lf]. 
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The BHs in QC-0 perform about half of an orbit prior to merging [131 E]; that is, QC-0 
looks more like a plunge/grazing collision. The intersection of the event horizon with 
the initial slice, however, has the topology of two separate spheres [TTj . 

The puncture BBH data of the Bowen-York type are defined by the bare masses 
7711,2 of the BHs, their coordinate locations C\$, assumed to be along the x-axis in the 
xy-plane, and their linear momentum P 12 , pointing along the y-axis. In the construction 
of the initial data, we vary m 2 = {1, 1.2, 1.3, 2, 4} mi while keeping mi fixed to its QC-0 
value of 0.453. We also keep fixed to the QC-0 values the puncture coordinate separation 
d = \C\ — C2I = 2.34 and the momentum parameters P\p — ±0.333, which means that 
the angular momentum value, J = 0.779, also remains unchanged. Note that the total 
ADM mass, Madm 5 of the configurations does change as do J and P when given in ADM 
mass units. Due to this setup, our initial data of unequal mass BBH systems do not obey 
the quasi-circular orbit condition of minimal binding energy [23]. The present work is 
aimed at investigating the effects of varying, in the QC-0 set-up, the bare mass parameter 
only. The motivation behind this choice was to study gravitational recoil starting from 
the IS CO plunge with the simplest parameter exploration. There is strong indication 
that the kick imparted to the BH that has resulted from the merger is dominated by the 
gravitational recoil during the plunge from ISCO |16j . We are currently extending the 
present study and investigating both unequal mass BBH mergers with initial separations 
outside of ISCO in quasi-circular orbit and plunge configurations with Post-Newtonian 
orbital parameters [16J. 

Table [1] summarizes the parameters in our simulations. We list the total mass 
M = Mi + M 2 , the mass ratio parameter 77 = M\.\ l>'M J where Mi j2 are the irreducible 
masses of the BHs computed from their individual apparent horizon (AH) well 
as angular momentum J and momentum parameter P in terms of the reduced mass, 
fj, = M\MifM. Table [1] also provides the time £ah in Max>m units when a common 
AH is first found. For QC-0 the orbital period of the equal mass case is estimated 
as t = 37.4 Madm- The drop in the time to merger t;ah in our results should not be 
interpreted as "unequal mass binaries merge faster". The effect is due to our approach 
in which the angular momentum of the configuration decreases as the bare mass ratio 
is decreased. Ideally one would like to compare initial data sets which are far separated 
and have the same orbital frequency. It is possible that for our cases with q < 0.55, 
a common event horizon already exists in the initial data slice; and, therefore, we are 
evolving a single, distorted BH. 

Methods. The evolutions were carried out using a code that solves the BSSN 
3+1 formulation of Einstein's equations [24[ [25j [26] . We use the "moving punctures" 
approach without excision [13], [13]. The code was produced by the Kranc code generation 
package [27] and uses the Cactus infrastructure. The simulations were performed using 
fourth order accurate centered finite differencing, except for the advection terms which 
were one-sided and second order accurate. The temporal updating is carried out with a 
three-step iterative Crank-Nicholson scheme with Courant factor of 0.25. Tests in the 
waveform for the equal-mass setup using resolutions h = {1/16,1/20,1/32} produced 



Unequal Mass Binary Black Hole Plunges and Gravitational Recoil 



4 



q 
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M 2 /M 
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P/fi 


J/(/iM) 


^AH 


1.00 


0.512 


0.512 


1.024 


0.250 


1.290 


2.948 


18.4 


0.85 


0.482 


0.565 


1.047 


0.248 


1.180 


2.637 


12.2 


0.78 


0.463 


0.597 


1.060 


0.246 


1.122 


2.476 


9.9 


0.55 


0.391 


0.712 


1.103 


0.229 


0.933 


1.979 


5.5 


0.32 


0.277 


0.872 


1.149 


0.183 


0.692 


1.409 


2.5 



Table 1. Initial data parameters and properties: The q = 1 case is the QC-0 data set 
in [22] • All models have puncture bare mass mi = 0.453, momenta P\ t 2 = ±0.333 and 
coordinate separation d = 2.34. The total mass M is given in Madm units as well as 
the coordinate time txa at which a common apparent horizon was found. 

convergence slightly below second order. Mesh refinement in the code is provided by 
Carpet [28], and tracking of AHs is handled by AHFinderDirect |29j . 

The gauge conditions used were modified versions of the 1+log lapse and T-driver 
shifts. Specifically, the lapse a was evolved using dtct = —2aK, where K is the trace 
of the extrinsic curvature. On the other hand, the shift vector was obtained from [14J: 
dtf3 l = FB % and dtB 1 = dtT 1 — fPdjY* — ^B l with £ a constant dissipative parameter 
and F = 3a/4, which guarantees that the asymptotic gauge speed associated with the 
longitudinal shift components is equal to the speed of light. The evolutions were started 
with (3 l and B % = 0. The advection term ftdjT 1 removes certain zero-speed modes of 
the system as analyzed in [30] . The parameter £ can be used to tune the rate of horizon 
expansion over the course of the evolution; large values lead to faster horizon growth. 
We have used values in the range 2 < £ < 5 and have not found any instabilities in this 
range. For the runs reported in this paper we used £ = 4. 

Our computational domain consisted of fixed 2:1 mesh refinements with 5 levels. 
The finest grid spanned —2 < x, y < 2 and < z < 2, with the coarsest —96 < x, y < 96 
and < z < 96 (we use bitant symmetry in the z direction). All the refinement levels 
except for the coarsest have the same number of grid points. We have run the unequal 
BBH models at two different resolutions, h = 1/16, 1/20. 

Results. Fig. [1] shows the irreducible mass of the AHs as a function of time, where 
Mi rr = a/ A/16it and A is the AH area. In all simulations with q > 0.32, we are 
able to track very accurately M irr for both of the individual merging BHs, as well as the 
resulting BH, over the entire course of the simulation. In the case of q = 0.32, we cannot 
track the smaller individual irreducible mass very accurately before merger. There is 
also a spurious drift in the common apparent horizon mass beyond the initial ADM 
mass content of the spacetime but the error only grows about 2% over t = 50 Madm- 
Notice that the smaller the value of q, the earlier the merger as indicated by tAu i n 
Table [TJ Post-merger, the individual AHs loose their meaning as apparent horizons 
since by definition the AH denotes the outermost marginally trapped surfaces. The 
individual surfaces remain marginally trapped. 

Fig. [2] shows snapshots of the horizons every 4.6 Madm before merger and at 
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Figure 1. The irreducible mass of the apparent horizon as a function of time for 
different mass ratios q. 
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Figure 2. Snapshots of the AH location in the xy-plane for the case q = 0.78. The 
larger BH is on the left moving toward the bottom. The snapshots are taken every 
4.6Madm prior to merger and at t = {40, 80, 105} Madm after merger. The first 
common AH is found at t = 9.9 Madm- The trajectories of the AH centroids are also 
shown. The common AH moves to the right and slightly upward after merger. 
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Figure 3. The dominant modes (£ = 2, m = 2 and I — 3, m — 3) of the real part of 
the Zerilli function tpi m against time for the different q ratios. The waveforms were 
extracted at r = 15. The £ — 2, m = 2 mode decreases in amplitude with decreasing q 
while the I = 3, m = 3 mode increases and then decreases again. 

t = {40, 80, 105} Madm after merger for the q = 0.78 case. The other cases are 
qualitatively similar. Notice how the initial common AH has an asymmetric peanut 
shape due to the unequal masses. Soon after it appears, the common AH becomes 
spherical, as the dynamical gauges drive the coordinates toward those of a single BH. 
At that moment, the common AH begins to drift slowly away from the origin. The last 
AH snapshot in Fig. [2] was taken at t = 105 Madm- The evident drift in the coordinate 
location of the common AH provides a hint that a kick is generated as a consequence 
of gravitational recoil. 

For waveform extraction, we compute Zerilli modes ip£ m using the Abrahams & Price 
convention |31j. Formally, the method assumes a spherically symmetric background. 
Simulations of QC-0 data [131 [14"] have produced rotating BHs with Kerr parameter 
a ~ 0.7. The main effect of using Zerilli extraction in spacetimes with significant 
angular momentum content will be a spurious amplitude offset in the waveform [321I33J. 
Fig. [3] displays the dominant modes (£ = 2, m = 2 and I = 3, m = 3) for the different 
q simulations. The extraction surface is located at r = 15 with the outer boundary at 
r = 96. Notice in Fig.[3]that the i = 2, m = 2 mode simply decreases in amplitude as the 
q ratio is decreased. This is mostly due to the initial data approaching that of a single 
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1.1 ±0.8 
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7.4 ±0.8 
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0.4 ±0.2 


0.011 ±0.006 


0.16 ±0.06 


2.6 ±0.6 


90 ±50 





32 


0.05 


0.001 


0.024 


0.4 
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Table 2. Estimates of radiated energy E la d (total), Eq (m = mode) and E2 (m = 2 
mode) as a % of the total initial ADM energy. Radiated angular momentum J ra d as a 
% of the initial angular momentum. Kick velocities V from gravitational recoil. The 
errors are the differences between the h = 1/16 and 1/20 resolution results multiplied 
by 2. 



distorted BH. The i = 3, m = 3 mode is smallest for the equal mass case (where it should 
be 0) and for the q = 0.32 case. Around t = 50Madm, the lower amplitude waveforms 
become affected by outer boundary effects. The strongly dominant i = 2, m = 2 mode 
remains accurate until the end of the simulation except for the q = 0.32 case where even 
the dominant mode is affected by the boundary by t = 60Madm- 

From the extracted modes, Table [2] summarizes estimates of the energy E T3 a and 
angular momentum J ra d emitted as % of the total energy and angular momentum as 
well as the recoil velocities V in km/s [31]. We compute the radiated energy from 
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and the recoil velocity from 

1 

V2r 



h. 



ihs 



E 



■0 



dp k 

~dt 



16tt 



Cm 



dh+ 
~dt 



+ 



dh x 
~dt 



lm 1 



n k dQ. 



Here ^^Y (9, 4>) denotes the spin-weight —2 tensor spherical harmonics. The kick 
velocity is obtained from MV l — J P l dt with M the total mass of the binary. The 
reported recoil velocity is \V\. 

As expected, the radiated energy and angular momentum are strongly dominated 
by the I = 2, m = 2 mode. The energy and angular momentum radiated for the q — 1 
case are in good agreement with previous work [TU [13]. The numbers reported in 
the table are computed from the h = 1/20 simulations. The errors are the difference 
between the quantities computed from the h = {1/16,1/20} resolution simulations 
multiplied by 2. The Richardson error estimate for our resolutions and 2nd order 
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Figure 4. Recoil velocity accumulated against time for the different models. Note 
that in comparison to far separated inspiral models ([MIH]) the inspiral contribution 
is missing. Also the q — 0.32 case shows an initial offset from zero which comes from 
the too close extraction radius. 

convergence is 16/9 and 2 is used as a conservative bound. We do not report errors 
for the q = 0.32 because the h — 1/16 simulation has insufficient resolution to compute 
the waveform. With increasing resolution the estimates of -E ra d and J ra d did increase, so 
assuming convergence, one would expect the continuum estimates for radiated energies 
and angular momenta to be slightly higher than those reported here. 

In order to compare the gravitational wave content of different simulations, Table [2] 
also presents the energy emitted in the m = and m = 2 modes, i.e. E and E 2 
respectively, as a % of the total initial ADM energy. The £ = 2, m = waveform in 
the q = 1 case changes quite strongly between the h = 1/20 and h = 1/16 simulations, 
resulting in large error bars on E Q and V (which should be zero for the equal mass 
case). Since the effect is not as strong for the £ = 2,m = 2 mode, the dominant mode 
for E and J ra d, the error bars are smaller. The £ = 2, m = 2 mode dominates for a pure 
inspiral, whereas the £ = 2, m = mode dominates for a plunge, so we would expect a 
larger E 2 /E ratio for inspirals than for plunges. It is reassuring that differences in the 
radiated energies of these modes exist; however, due to the construction of the initial 
data sequence, these differences do not reflect a change of q only. It will be interesting 
to see how this changes with further separated, truly inspiraling configurations starting 
from outside ISCO. 

The recoil velocities, V, reported in Table [1] are consistent with those obtained 
from head-on collisions [37], mixed numerical-perturbative inspiral mergers [35] and 
recent fully relativistic results at larger separations [36, 18J. For reference, kicks from 
gravitational recoil of relevance to galactic BH merger scenarios [9] have been recently 
estimated to second post-Newtonian order [IB] . The kicks were found to be dominated 
by the plunge phase and could reach speeds larger than 100 km/s for 0.1 < n < 0.24 
or 1/8 < q < 2/3. Recent work using the effective one-body approach gives lower kick 
velocities of at most 74 km/s at rj = 0.2 or q = 0.38 (TTJ. 
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reduced mass parameter r| 

Figure 5. A comparison of recoil velocity estimates as a function of 77. The 
2PN estimates of Blanchet et al. [TB] are denoted by the solid line and those from 
Gopakumar-Damour [T7] by a dashed line. Campanelli [35] and Baker et al. [35] kicks 
are labeled by a circle and a box, respectively, and our recoil velocities are labeled with 
diamonds. Note that the smaller r\ cases present lower bounds due to the initial data 
configurations studied. 



In Fig. H] we show the recoil velocity accumulated as a function of time. The close 
detector radius manifests itself in the non-zero initial offset of the recoil velocity of 
the q = 0.32 model. The inspiral contribution is notably absent from these models 
due to their close separation (compare Refs. [361 [18]). Also note that the initial data 
contains an initial feature that we cannot remove as it is already overlayed with the 
merger signal. In further separation models as in Ref. [18] the initial feature shows up 
in the weaker inspiral part and can therefore easily be removed. Fig. [5] shows our recoil 
velocities as a function of the reduced mass parameter 77, including kick estimates found 
by others pH QH EH [36] . Our q = 0.85, n = 0.248 and q = 0.78, n = 0.246 cases are 
comparable to the estimates from Blanchet et al. [16]. The reason for the agreement is 
that these cases are closer to the QC-0 equal mass case. Thus the initial data setup is 
not too far from ISCO and quasi-circularity. The case q = 0.55, rj = 0.229 yields the 
largest recoil velocity, ~ 82 km/s, also consistent with Blanchet et al. PS] who find a 
peak in the recoil around q = 0.4, n = 0.204. Recently Baker et al. [36] reported a recoil 
velocity of 105 km/s for q = 0.67, n = 0.24, a value larger than our kicks. The difference 
is mostly due to the larger initial separation of their BHs. Our smallest kick is obtained 
in the q = 0.32, n = 0.183 case. This is likely a significant underestimation since we are 
probably dealing with a single distorted BH. 

Conclusions. We have shown results from a series of unequal mass BBH inspiral 
simulations. The kick velocities spanned a range of 25-82 km/s. Our kick velocity 
estimates should be understood as a lower bound because the initial separations are 
smaller than those normally used (for example in Ref. [16]) for the plunge phase. The 
observed drop in the time to merger £ah in our results comes from our approach in 
which the angular momentum of the configuration decreases as the bare mass ratio is 
decreased. This approach introduces significant uncertainty in the recoil extracted and 
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in particular it is possible that for our cases with q < 0.55, a common event horizon 
already exists in the initial data slice; and, therefore, we are evolving a single, distorted 
BH. Ideally one would like to compare initial data sets which are far separated and 
have the same orbital frequency. We would also like to point out that these close 
configurations contain significant eccentricities which are known to influence the recoil 
velocities. For small eccentricities e < 0.1 Ref. [2] shows Vr oc (1 + e). More on equal 
mass inspirals and eccentricities can be found in Ref. [381 139] . 

In addition to the gravitational recoil, we estimated the radiated energy and angular 
momentum, paying particular attention at the energy radiated independently in the 
m = 0, 2 modes. We also monitored the irreducible masses during the simulations, 
providing a good indicator that the near-zone physics was accurately evolved. A more 
detailed comparison of the dependence of the waveforms on the mass ratio from inspirals 
starting outside the ISCO is currently under investigation. A strong dependence of the 
waveform on the binary parameters would facilitate their estimation, but at the same 
time it would hinder initial detection efforts as many templates would be required. If 
waveforms, on the other hand, do not vary significantly with mass ratios, the search 
effort could be easier at the expense of accurate parameter estimation. 
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